
defsysv,'!pc', EXIST=exist
IF ~exist THEN omni_physical_constants

t = read_mrt('./ancillary/RomanDuval_2010_Table1.mrt')
ancil = omni_read_conffile('./conffiles/ancillary.conf')

help,t
help,t,/str

vvir = sqrt(2.d * !pc.G * t.mass*!pc.msun / (t.r * !pc.pc) ) / 1.d3 ; km/s

myps,'./grs_paper/grs_dv_vvir.eps',xsize=10,ysize=10
multiplot_xm,[2,2],mpcharsize=1.0,xgap=0.035,/hardleft,ygap=0.035,/doxaxis

cgPlot,vvir,t.delv,charsize=1.0,psym=16,symsize=0.4,xtit='Virial  [km/s]',$
       ytit='Measured  [km/s]'

one2one,color='cyan',linestyle=3

linfit_plot,vvir,t.delv,/legend,/top,/left,color='red'

multiplot,/doyaxis,/doxaxis

alpha = (t.delv / vvir) * (t.delv / vvir)

cgPlot,t.alpha,alpha,charsize=1.0,psym=16,symsize=0.4,xtickformat='exponent10',$
       xtit='Virial Parameter',ytit='(DELV / VVIR)!u2!n',/xlog,/ylog,$
       ytickformat='exponent10'
one2one,color='cyan',linestyle=4,/log


print,m4_stat(alpha/t.alpha)

multiplot,/doyaxis,/doxaxis

bin = 0.25
plothist,vvir,bin=bin,xtit='Velocity Dispersion  [km/s]',$
         ytit='N per '+string(bin,format="(F0.2)")+' km/s bin',$
         charsize=1.0,yr=[0,100]
plothist,t.delv,bin=bin,/over,color='cyan'

al_legend,/top,/right,box=0,linestyle=0,color=['cyan','black'],linsize=0.6,$
          ['Measured','Virial']


multiplot,/doyaxis,/doxaxis

plothist,alog10(t.alpha),bin=0.1,xtit='Log!d10!n Virial Parameter',$
         charsize=1.0,xst=4
cgAxis,xaxis=0,xr=10.^!x.crange,/xst,charsize=1.0,xtit='Virial Parameter',$
       /xlog,/save,xtickformat='exponent10'
cgAxis,xaxis=1,xtickformat='blank_axis',/xst,/xlog

myps,/done,/mp


val = sqrt(!pc.G * 3.0d4*!pc.msun / (1.d * !pc.pc) ) / 1.d3 ; km/s
print,val,val / sqrt(3.d)


print,1.d*!pc.pc*(7.0d3*7.0d3*3.d)/!pc.G/!pc.msun

;;======================================================================


myps,'./grs_paper/grs_phys_cdf.eps',xsize=9,ysize=3.5

multiplot_xm,mpcharsize=2,[2,1],/hardleft

rad = t.r
print,'RAD: ',n_elements(rad)
rcdf = ccdf(rad,/cdf)

cgPlot,rad,rcdf,charsize=1.0,xtit='Physical Radius  [pc]',$
       ytit='Cumulative Distribution (P<)',xr=[0,34.999]

rval = interpol(rad,rcdf,0.9)
vline,rval,color='cyan'
print,rval
vline,median(rad),color='red'
print,median(rad)
rval2 = interpol(rad,rcdf,0.1)
vline,rval2,color='green'
print,rval2
al_legend,box=0,/bottom,/right,charsize=0.8,$
          ['10','50','90']+'%  have R < '+$
          string([rval2,median(rad),rval],format="(F0.1)")+' pc'


multiplot

vel = t.delv

vcdf = ccdf(vel,/cdf)

cgPlot,vel,vcdf,charsize=1.0,xtit='FWHM Velocity Dispersion  [km s!u-1!n]'
       
vval = interpol(vel,vcdf,0.9)
vline,vval,color='cyan'
print,vval
vline,median(vel),color='red'
print,median(vel)
vval2 = interpol(vel,vcdf,0.1)
vline,vval2,color='green'
print,vval2
al_legend,box=0,/bottom,/right,charsize=0.8,$
          ['10','50','90']+'% have '+cgSymbol('Delta')+'V < '+$
          string(vval2,median(vel),vval,format="(F0.1)")+' km s!u-1!n'
vline,median(vel),color='red'



myps,/done,/mp


;;========================================================
;;========================================================
;;========================================================
;;========================================================
;;========================================================
;; 2-D Plots
myps,'./grs_paper/grs_phys_rv.eps',xsize=13
multiplot_xm,[3,1],mpcharsize=1.0,xgap=0.001,/xmargin

;; cgPlot,t.r,t.delv,psym=16,symsize=0.4,xtit='Physical Radius  [pc]',$
;;        ytit='FWHM Velocity Dispersion  [km s!u-1!n]',charsize=1.0,xr=[0,35],$
;;        xtickinterval=10
;; vline,16.5,color='cyan'
;; vline,/h,5.1,color='cyan'
;; ind = where(t.r LE 16.5 AND t.delv LE 5.1, nle)
;; print,'N within 90% BOTH: ',nle,double(nle)/n_elements(t)

;; vline,7.4,color='brown'
;; vline,/h,3.1,color='brown'
;; ind = where(t.r LE 7.4 AND t.delv LE 3.1, nle)
;; print,'N within 50% BOTH: ',nle,double(nle)/n_elements(t)

;; pcdf = dblarr(n_elements(t))
;; FOR ii=0,n_elements(t)-1 DO BEGIN
;;    ind = where(t.r LE t[ii].r AND t.delv LE t[ii].delv, nle)
;;    pcdf[ii] = double(nle) / n_elements(t)
   
;; ENDFOR

;; levels = findgen(11)/10.
;; print,levels
;; cgContour,pcdf,t.r,t.delv,/irregular,levels=levels,color='red',/over,$
;;           resolution=[51,51]


;; multiplot

nr = 200
nv = 200

rrange = [0.5,30.5]
vrange = [1.9,9.3]

vmax = vrange[1] - vrange[0]
rmax = rrange[1] - rrange[0]

rarr = findgen(nr)*rmax/nr + rrange[0]
varr = findgen(nv)*vmax/nv + vrange[0]
cdf2d = dblarr(nr,nv)
rr = cmreplicate(rarr,nv)
vv = transpose(cmreplicate(varr,nr))

FOR ii=0,nr-1 DO BEGIN
   FOR jj=0,nv-1 DO BEGIN
      
      ind = where(t.r LE rarr[ii] AND t.delv LE varr[jj], nle)
      cdf2d[ii,jj] = double(nle) / n_elements(t)
   ENDFOR
ENDFOR

cgLoadct,0
cgPlot,/nodata,[0],[0],title='GRS Cloud Bivariate CDF',charsize=1.0,xst=4,yst=4
cgLoadct,75
plotimage,cdf2d,range=double([0,1]),xst=4,yst=4,min_dpi=100


levels = [0.1,0.25,0.5,0.75,0.9,1.0]

cgContour,cdf2d,levels=levels,color='black',/over,charsize=1.3,$
          c_annot=string(levels,format="(F0.2)")


;; arad = sqrt(rr*rr/(rmax*rmax)+vv*vv/(vmax*vmax))
;; cgContour,arad,/over,color='violet'

;; ind = where( abs(cdf2d - 0.75) LE 0.01 AND arad LE 0.65, ni)
;; print,'-----'
;; print,m4_stat(rr[ind],vv[ind],cdf2d[ind],arad[ind])
;; print,'-----'


cgAxis,xaxis=0,xr=rrange,/save,xtit='Physical Radius  [pc]',charsize=1.0,$
       /xst,xtickinterval=10
cgAxis,xaxis=1,/xst,xtickformat='blank_axis'
cgAxis,yaxis=0,yr=vrange,/save,ytit='FWHM Velocity Dispersion  [km s!u-1!n]',$
       charsize=1.0,/yst,yminor=4
cgAxis,yaxis=1,/yst,ytickformat='blank_axis',yminor=4

cgPlots,psym=14,color='brown',symsize=2,  16.5, 5.1
cgPlots,psym=4,color='black',symsize=2,   16.5, 5.1
cgPlots,psym=14,color='bisque',symsize=2,  7.4, 3.1
cgPlots,psym=4,color='black',symsize=2,    7.4, 3.1
cgPlots,psym=14,color='blu3',symsize=2,   14.7, 4.75
cgPlots,psym=4,color='black',symsize=2,   14.7, 4.75
cgPlots,psym=14,color='pink',symsize=2,   10.5, 3.6
cgPlots,psym=4,color='black',symsize=2,   10.5, 3.6
cgPlots,psym=14,color='grn4',symsize=2,    9.0, 4.1
cgPlots,psym=4,color='black',symsize=2,    9.0, 4.1


multiplot;,/doyaxis

;;&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

ndr = fix(floor(4*5*1.5))
ndv = fix(floor(5*5*1.5))

dr = (dindgen(ndr)*1.d + 1.d)*1.d
dv = dindgen(ndv)*0.2d + 2.d

restore,'./grs_paper/assoc_volume.sav',/ver
contcol='black'
print,m4_stat(nhr)
cgLoadct,0
cgPlot,/nodata,[0],[0],title='BGPS Objs w/ '+cgSymbol('geq')+'1 HRDS HII Region',charsize=1.0,xst=4,yst=4
cgLoadct,75
pr = set_plot_range(nhr)
print,'N HRDS',pr
plotimage,nhr,xst=4,yst=4,range=pr,min_dpi=100
cgLoadct,0

ddr = arrstep(dr)
ddv = arrstep(dv)
print,'X Axis: ',[min(dr),max(dr)]+[-1,1]*ddr/2.
print,'Y Axis: ',[min(dv),max(dv)]+[-1,1]*ddv/2.
cgAxis,xaxis=0,xr=[min(dr),max(dr)]+[-1,1]*ddr/2.,/xst,charsize=1.0,/save,$
       xtit='Physical Radius  [pc]'
cgAxis,xaxis=1,/xst,xtickformat='blank_axis'
cgAxis,yaxis=0,yr=[min(dv),max(dv)]+[-1,1]*ddv/2.,/yst,charsize=1.0,/save,$
       ytickformat='blank_axis',yminor=4
cgAxis,yaxis=1,/yst,ytickformat='blank_axis',yminor=4

cgContour,nhr,dr,dv,levels=findgen(21)*200.,/over,color=contcol,charsize=1.3

cgPlots,psym=14,color='brown',symsize=2,  16.5, 5.1
cgPlots,psym=4,color='black',symsize=2,   16.5, 5.1
cgPlots,psym=14,color='bisque',symsize=2,  7.4, 3.1
cgPlots,psym=4,color='black',symsize=2,    7.4, 3.1
cgPlots,psym=14,color='blu3',symsize=2,   14.7, 4.75
cgPlots,psym=4,color='black',symsize=2,   14.7, 4.75
cgPlots,psym=14,color='pink',symsize=2,   10.5, 3.6
cgPlots,psym=4,color='black',symsize=2,   10.5, 3.6
cgPlots,psym=14,color='grn4',symsize=2,    9.0, 4.1
cgPlots,psym=4,color='black',symsize=2,    9.0, 4.1

multiplot,/doxaxis,/doyaxis


cgLoadct,0
cgPlot,/nodata,[0],[0],title='HRDS vs. EMAF + H!d2!n as Prior DPDF',$
       charsize=1.0,xst=4,yst=4

fma = median(fma,3)
print,m4_stat(fma)
cgLoadct,75
pr = set_plot_range(fma)
print,'Matching: ',pr
plotimage,fma,xst=4,yst=4,range=pr
cgLoadct,0
ddr = mean(dr[1:*] - dr[0:*])
ddv = mean(dv[1:*] - dv[0:*])

cgAxis,xaxis=0,xr=[min(dr),max(dr)]+[-1,1]*ddr/2.,/xst,charsize=1.0,/save,$
       xtit='Physical Radius  [pc]'
cgAxis,xaxis=1,/xst,xtickformat='blank_axis'
cgAxis,yaxis=0,yr=[min(dv),max(dv)]+[-1,1]*ddv/2.,/yst,charsize=1.0,/save,$
       ytickformat='blank_axis',yminor=4
cgAxis,yaxis=1,/yst,ytickformat='blank_axis',yminor=4


cgContour,nehr,dr,dv,levels=findgen(21)*50,/over,color=contcol,charsize=1.3
cgContour,nhr,dr,dv,levels=findgen(21)*200.,/over,color='brown',charsize=1.3
cgContour,fma,dr,dv,levels=[.8,0.85,0.9,0.95],/over,color='black',charsize=1.3
cgContour,cdf2d,rarr,varr,levels=levels,color='cyan',/over,charsize=1.3,$
          c_annot=string(levels,format="(F0.2)")


cgPlots,psym=14,color='brown',symsize=2,  16.5, 5.1
cgPlots,psym=4,color='black',symsize=2,   16.5, 5.1
cgPlots,psym=14,color='bisque',symsize=2,  7.4, 3.1
cgPlots,psym=4,color='black',symsize=2,    7.4, 3.1
cgPlots,psym=14,color='blu3',symsize=2,   14.7, 4.75
cgPlots,psym=4,color='black',symsize=2,   14.7, 4.75
cgPlots,psym=14,color='pink',symsize=2,   10.5, 3.6
cgPlots,psym=4,color='black',symsize=2,   10.5, 3.6
cgPlots,psym=14,color='grn4',symsize=2,    9.0, 4.1
cgPlots,psym=4,color='black',symsize=2,    9.0, 4.1

cgLoadct,75
cgColorbar,range=pr,position=[0.94,!y.window[0],0.96,!y.window[1]],$
           title='Distance Agreement Fraction',charsize=0.8,/vert,/right
cgLoadct,0





myps,/done,/mp

;;============================================================
;; Try the CCDF
myps,'./grs_paper/grs_phys_ccdf.eps',xsize=9,ysize=3.5

multiplot_xm,mpcharsize=1.0,[2,1];,/hardleft

rad = t.r

rcdf = ccdf(rad)

cgPlot,rad,rcdf,charsize=1.0,xtit='Physical Radius  [pc]',$
       ytit='CCDF',xr=[0,34.999]

rval = interpol(rad,rcdf,0.1)
vline,rval,color='cyan'
print,rval
vline,median(rad),color='red'
print,median(rad)
rval2 = interpol(rad,rcdf,0.9)
vline,rval2,color='green'
print,rval2
al_legend,box=0,/bottom,/right,charsize=0.8,$
          ['90','50','10']+'%  have R > '+$
          string([rval2,median(rad),rval],format="(F0.1)")+' pc'


multiplot

vel = t.delv

vcdf = ccdf(vel)

cgPlot,vel,vcdf,charsize=1.0,xtit='FWHM Velocity Dispersion  [km s!u-1!n]'
       
vval = interpol(vel,vcdf,0.1)
vline,vval,color='cyan'
print,vval
vline,median(vel),color='red'
print,median(vel)
vval2 = interpol(vel,vcdf,0.9)
vline,vval2,color='green'
print,vval2
al_legend,box=0,/bottom,/right,charsize=0.8,$
          ['90','50','10']+'% have '+cgSymbol('Delta')+'V > '+$
          string(vval2,median(vel),vval,format="(F0.1)")+' km s!u-1!n'
vline,median(vel),color='red'



myps,/done,/mp




;;&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
;;&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
;;&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
;; Make single-panel plot for the Association Volume Plot!

myps,'./grs_paper/assoc_volume_optimize.eps',xsize=3.1,ysize=2.8,thick=1.5,$
     ct=34
!x.omargin = [-0.75,+2.5]
!y.omargin = [-0.5,-1.5]
multiplot,[1,1]

csf = 0.55

cgLoadct,34;75

matlevel = [0.8,0.85,0.90]

pr = set_plot_range(nhr)
plotimage,nhr,range=pr,xst=4,yst=4,min_dpi=300


;;************************************************
;; Bivariate CDF Overplot!
nr = 200
nv = 200
rrange = [0.5,30.5]
vrange = [1.9,9.3]
vmax = vrange[1] - vrange[0]
rmax = rrange[1] - rrange[0]
rarr = findgen(nr)*rmax/nr + rrange[0]
varr = findgen(nv)*vmax/nv + vrange[0]

dx = arrstep(rarr)
dy = arrstep(varr)
cgAxis,xaxis=0,/xst,xr=rrange + dx/2.*[-1,1],/save,$
       xtit='Physical Radius  [pc]',charsize=1.0*csf
cgAxis,xaxis=1,/xst,xtickformat='blank_axis'
cgAxis,yaxis=0,/yst,yr=vrange + dy/2.*[-1,1],/save,$
       ytit='FWHM Velocity Dispersion  [km s!u-1!n]',charsize=1.0*csf
cgAxis,yaxis=1,/yst,ytickformat='blank_axis'

;; Number Contours
cgContour,nhr,dr,dv,levels=[200,400,600,800,1000,1200,1400],$
          /over,color='blk7',c_charsize=1.0*csf,c_linestyle=5

levels = [0.1,0.25,0.5,0.75,0.9,1.0]
cdf2d = median(cdf2d,5)
cgContour,cdf2d,rarr,varr,levels=levels,color='black',/over,c_charsize=1.0*csf,$
          c_annot=string(levels,format="(F0.2)")
;;************************************************


;; Matching percentage
;; cgContour,fma,dr,dv,levels=matlevel,/over,c_charsize=1.0*csf,$
;;           c_color='deep pink',c_annot=string(matlevel*100.,format="(' ',I0,'% ')")
fma = median(fma,3)
cgContour,fma,dr,dv,levels=matlevel,/over,c_charsize=1.0*csf,$
          c_color='wt1',c_annot=string(matlevel*100.,format="(' ',I0,'% ')")

help,fma,dr,dv



cgPlots,psym=14,color='brown',symsize=1,  16.5, 5.1
cgPlots,psym=4,color='black',symsize=1,   16.5, 5.1
cgPlots,psym=14,color='blk5',symsize=1,  7.4, 3.1
cgPlots,psym=4,color='black',symsize=1,    7.4, 3.1
;; cgPlots,psym=14,color='blu3',symsize=2,   14.7, 4.75
;; cgPlots,psym=4,color='black',symsize=2,   14.7, 4.75
;; cgPlots,psym=14,color='pink',symsize=2,   10.5, 3.6
;; cgPlots,psym=4,color='black',symsize=2,   10.5, 3.6
;; cgPlots,psym=14,color='blk5',symsize=1,    9.0, 4.1
;; cgPlots,psym=4,color='black',symsize=1,    9.0, 4.1
cgPlots,psym=14,color='white',symsize=1,ancil.ppv_dis,ancil.ppv_dv
cgPlots,psym=4,color='black',symsize=1,ancil.ppv_dis,ancil.ppv_dv

cgColorbar,/vert,/right,position=[0.87,!y.window[0],0.90,!y.window[1]],$
           range=pr,charsize=0.9*csf,$
           title='BGPS Objs w/ '+cgSymbol('geq')+'1 HRDS HII Region'

myps,/done

END
